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A first principles embedded cluster approach is used to calculate O chemical shielding tensors, 
a, in prototypical transition metal oxide ABO3 perovskite crystals. Our principal findings are 1) a 
large anisotropy of a between deshielded — Oy and shielded cr^ components (2 along the Ti-O 
bond); 2) a nearly linear variation, across all the systems studied, of the isotropic aiso and uniaxial 
CTax components, as a function of the B-O-B bond asymmetry. We show that the anisotropy and 
linear variation arise from large paramagnetic contributions to Ox and Oy due to virtual transitions 
between 0(2p) and unoccupied B(nd) states. The calculated isotropic Siso and uniaxial Jax chemical 
shifts are in good agreement with recent BaTiOa and SrTiOa single crystal ^'^O NMR measurements. 
In PbTiOs and PbZrOa, calculated 5iso are also in good agreement with NMR powder spectrum 
measurements. In PbZrOa, (5iso calculations of the five chemically distinct sites indicate a correction 
of the experimental assignments. The strong dependence of a on covalent 0(2p)-B(nd) interac- 
tions seen in our calculations indicates that ^^O NMR spectroscopy, coupled with first principles 
calculations, can be an especially useful tool to study the local structure in complex perovskite 
alloys. 
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I. INTRODUCTION 



Due to their reversible polarization and strong electro- 
mechanical coupling, ferroelectric perovskite oxides are 
key components in many electronic and mechanical de- 
vices such as sensors, actuators, and random access 
mcmoryji The largest piezoelectric response occurs when 
competing structural instabilities are present^ so that 
the local structure depends sensitively on strain, exter- 
nal fields, and chemical composition, as for example, 
near the morphotropic phase boundary near x = 1/2 
in Pb(Zri_a;Tia;)03 (PZT). Nuclear magnetic resonance 
(NMR) can be an important probe of local structure, and 
NMR has increasingly been used to study complex ferro- 
electric alloysi^i^ First-principles calculations of NMR 
properites are likely to play an important role in helping 
to interpret complex measured spectra. 

A material's characteristic NMR spectra is determined 
by the coupling of the nuclear magnetic dipole and elec- 
tric quadrupole moments with the local magnetic field 
and with the electric field gradients (EFG), respectivelyi^ 
The chemical shielding tensor, ct, determines the local 
magnetic field at a nucleus 



(PBC) 



14.15.16 



and implementations have been limited to 



B = (1 - (7)Bext, 



(1) 



where the induced field Bind = ^o'Boxt arises from elec- 
tronic screening currents. Theoretical determination of 
the chemical shielding tensor, &, is more subtle than 
the EFGs, since the latter depend only on the ground 
state charge density. Linear response methods have 
long been used to calculate a for isolated molecules and 
clusters jii^iiSiiiSiiiii^ii^ In extended systems complications 
arise due to the use of periodic boundary conditions 



planewave-based methods. Recently, an alternative to 
the linear response method has been proposed based on 
calculations of the orbital magnetization] If this 
approach proves successful, calculations of a could be 
implemented in other standard band structure codes, in- 
cluding all-electron methods 

Calculations of & in crystals and extended systems 
have largely used finite-size clusters, due to the lack of 
generally available PBC computer codes with NMR func- 
tionality. Accurate results can be obtained with this ap- 
proach if the screening currents are sufficiently localized 
near the target nucleus, and if its local environment is 
adequately modeled by the cluster. Embedded cluster 
techniques have been successfully used in crystals and 
macromolecules to obtain nuclear quadrupole resonance 
spectrum, ligand to metal charge transfer excitations, 
photoemission, electric field gradients, and hyperfine cou- 
pling in high Tc superconductorsi ^^i^^'^'^i^^ 

Such calculations have generally used standard all- 
electron quantum chemistry methods, which are widely 
available in many quantum chemistry computer pro- 
grams, such as GAUSSIANj^iS^ These methods are very 
mature and can calculate a with a range of approxi- 
mations, using well tested Gaussian type orbital (GTO) 
basis setsi^ In increasing order of computational cost, 
these methods range from mean-field type [density func- 
tional theory (DFT), Hartree Fock (HF), and hybrid 
methods] to correlated approximations [such as second- 
order MoUer-Plesset (MP2) perturbation and coupled 
cluster (CC) methods] i^ '^^i^^ The present cluster calcula- 
tions exploit this flexibility to calculate chemical shield- 
ings, comparing DFT/GGA, restricted HF (RHF), and 
hybrid-DFT/B3LYP calculations. Currently, only LDA 
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and GGA DFT approximations are available for NMR 
calculations in PBC methods. 

For covalently bonded systems, cluster calculations 
typically use hydrogen atoms to terminate dangling 
bonds at the cluster surfaced Instead, we use point 
charge embedding2£ to model the long range Coulomb 
interactions in the ABO3 materials studied here. Ad- 
ditional techniques are also used to better handle the 
polarizable character of the large O^^ anion and also to 
control surface depolarization effects. 

In this paper, we show that the embedded cluster ap- 
proach can yield accurate oxygen a for transition metal 
perovskites. The focus is on four prototypical materials: 
BaTiOa, SrTiOg, PbTiOa and PbZrOa- Strong cova- 
lency between the O and transition metal B atoms (and 
between the A and O atoms in Pb based materials) del- 
icately balance ionic electrostatic interactions in these 
materials and related alloys, resulting in a wide variety 
of interesting and technologically important properties. 
The classic ferroelectric material BaTiOs exhibits three 
reversible temperature dependent ferroelectric phases. 
SrTiOa is a key constituent in superlattice structures 
with novel material propertiesi^ PbTiOa is the proto- 
typical Pb based ferroelectric. The strong covalency of 
Pb on the A-site is responsible for the large c/a = 1.065 
ratio in PbTi03, and it plays a similar role in large elec- 
tromechanical response of relaxor ferroelectrics such as 
PMN.— PbZrOs adopts a complicated non-polar antifer- 
rodistortive structure with five chemically inequivalent O 
sites, which our calculations reproduce. PbZrOa is also 
of the end-point compounds, together with PbTiOa, of 
the PZT solid solution series. 

The paper is organized as follows. The theoretical ap- 
proach is described in Sec. [TTl Calculated results for 
ABO3 systems are presented in Sec. IIIIl and an anal- 
ysis of the results is presented in Sec. |TVl Finally, we 
summarize and conclude in Sec. |Vl 



II. THEORETICAL APPROACH 

After briefly reviewing g conventions, the construction 
of embedded clusters for ionic materials is described. Fi- 
nally, we describe the quantum chemistry methods and 
GTO basis sets used in the calculations. 



A. Chemical shielding tensor 

The chemical shielding tensor g [Eq. ([T])] is a mixed 
second derivative of the ground state energy with re- 
spect to the applied magnetic field and nuclear magnetic 
moment ii As such, it is an asymmetric second rank tensor 
with nine independent components in general, although 
its symmetry can be higher, depending on the site sym- 
metry of the target nucleusj^ The anti-symmetric part of 
a contributes negligibly to the NMR resonance frequency 
shift, since it enters only in second order ^^^i^ although it 



can contribute to relaxation.^ The symmetric part can 
always be diagonalized, and the NMR frequency is deter- 
mined by the following combinations of its principal axis 
components'^ 

Ciso I (o-x + + 0-2) = jTr G 
fax ~ \ (2(Tz - CTx - O-y) — \{<^z— CTiso) (2) 
f^aniso 2 ^x) : 

where (Xiso, Cax, and (Janiso are the isotropic, uniaxial, and 
anisotropic components, respectively. In this paper, g^ is 
chosen to correspond to the principal axis that is mostly 
nearly along the B-O-B bond direction. 

Positive values of g are conventionally denoted as 
shielding the external field, while negative elements are 
referred to as deshielding [see Eq. Measurements of 
G are usually reported with with respect to a reference 
compound, where the chemical shift tensor is defined as^^ 

S=-{g-gI^^). (3) 



B. Embedded clusters 

An embedded cluster model is depicted in Fig. [1] Its 
core consists of real, or "quantum" (QM) atoms, with 
the target O atom at its center. The total number of 
electrons is determined by satisfying the nominal formal 
valence of all the QM atoms. The A4B2O15 QM cluster in 
the figure thus has a net charge of -14. The QM cluster is 
embedded in the classical potential due to point charges, 
pseudopotentials, and, in some cases, an external electric 
field to cancel surface depolarization effects, as described 
in following subsections. The atomic site designations 
nn and nnn, used below, denote nearest-neighbor and 
next-nearest-neighbor sites, respectively, and are based 
on the ideal structures. Thus, for example, the Ti atom 
in tetragonal P4mm PbTiOs would be regarded as hav- 
ing six nn O atoms, despite the distortion of the TiOg 
octahedron. In the actual calculations, the tetragonal 
distortion results in two chemically inequivalent Ooq and 
Oax atoms (Oax lying on the polar axis), which require 
two separate embedded cluster calculations. 



1. Quantum cluster 

In all the QM clusters considered in this paper, the tar- 
get oxygen atom is fully coordinated with QM atoms lo- 
cated at its nn and nnn sites. Secondly, the target atom's 
nn QM atoms are themselves fully coordinated with nn 
QM atoms. Finally additional QM atoms are added, as 
required by ideal perovskite symmetry This procedure 
results in a 21 QM-atom cluster: (A4B20i,5)^"'~ , depicted 
in Fig. [H where A = Sr, Ba, or Pb; B = Ti or Zr. This 
includes two corner-shared BOg octahedra, centered on 
the targeted O atom (11 O and 2 B atoms), the target 
atom's 4 nnn A cations, and 4 additional O atoms, which 
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are also one lattice constant away from the center atom. 
Tests with larger QM clusters were used to estimate con- 
vergence of the chemical shielding tensor with respect to 
cluster size, as discussed in Sec. Ill CI 



2. Madelung potential: point charges 

Next, the QM cluster is embedded in the crystal envi- 
ronment by surrounding it with a large array of point 
charges. The purpose of the point charges is to bet- 
ter simulate the crystal environment by generating the 
correct crystalline electrostatic Madelung potential in 
the QM region. The correct Madelung potential also 
plays a key role in stabilizing the ion, as has been 
shown in Gordon-Kim^ type modelsj^i^ where the in- 
ternal energy of ionic systems is determined from the 
energy of overlapping ionic charge densities. The finite 
point charge distribution is determined using the EWALD 
program^ as follows. In a first step, EWALD calculates 
the Madelung potential with the Ewald method for PBC, 
using nominal ionic values {e.g., Qi ~ —2 and Qi = +2 
for and Pb^"*", respectively) for the atoms placed at 
crystallographic positions of the targeted system. In a 
second step, EWALD retains the nearest O(IO^) Qi cen- 
tered on the target atom, adjusting the values of the 
outermost Qi to reproduce the Madelung electrostatic 
potential on and in the vicinity of the QM atoms. In 
this second step, the nearest ~ 500 — 750 Qi are fixed at 
their nominal values, and, in addition, the net monopole 
and dipole moments of the point charge distribution are 
constrained to vanish. 



3. Boundary compensation via empty pseudopotentials 

To accelerate the convergence of a with respect to the 
size of the QM cluster, it is advantageous to control the 
artificial polarization of boundary 0(2p) states. This 
arises from the strongly attractive electrostatic poten- 
tial of neighboring cation point charges for QM oxygen 
atoms on the periphery of the cluster. To alleviate this, 
the nn and nnn cation point charges of boundary O atoms 
are replaced by "empty" pseudopotentials (ePSP))^^ as 
illustrated in Fig. [TJ The ePSP is defined as follows: 
i) it is a large core pseudopotential for the most loosely 
bound valence electrons of the corresponding cation {e.g., 
a Ti'^^ PSP); ii) there are no GTO basis functions asso- 
ciated with this site. The resulting modified cation clas- 
sical potential simulates the Pauli cation core repulsion, 
thereby reducing the artificial polarization of boundary 
0(2p) states4S In the embedded 21 atom QM cluster 
discussed above, for example, 26 boundary point charges 
are replaced by cPSPs (labeled as A* or B*), yielding a 
(A4B20i5)^"'~ - A*gB*Q cluster, which is surrounded by 
the remaining ~ 0(10^) point charges. 



4. Cancellation of electric depolarization fields 

Due to the spontaneous electric polarization found in 
some materials, such as tetragonal PbTiOs, a macro- 
scopic depolarizing electric field is present in finite sam- 
ples, in the absence of surface compensating charges. 
Calculations using PBC are implicitly done in zero total 
macroscopic electric field, which automatically excludes 
surface depolarization effects. In the present finite size 
clusters, a depolarizing electric field can arise from a pos- 
sible net dipole moment, due to polarization of the quan- 
tum mechanical charge density, i.e., from the wave func- 
tions of the QM cluster. The net dipole moment due to 
the point charges and ePSPs is zero by construction, as 
discussed above. The resulting depolarizing electric field 
is removed in the calculation, by applying an external 
electric field in the opposite direction. The magnitude of 
the external field is chosen so that the force on the cen- 
tral target atom matches that of an all-electron PBC lin- 
earized augmented planewave (LAPW)^ calculation. In 
normal equilibrium conditions, using experimentally de- 
termined structures, the LAPW forces are usually small, 
since the theoretical structure is usually close to that of 
experiment. 



C. Calculational details 

The embedded cluster calculations were performed 
with the GAUSSIAN computational package^^ The 
chemical shielding tensor was determined using the 
continuous set of gauge transformations (CSGT) 
methodi^i^ The gauge-independent atomic orbital 
(GIAO) method^ was also used for comparison in some 
cases. The gauge origin in GIAO is at the target atom, 
yielding a useful decomposition into diamagnetic and 
paramagnetic components, as discussed further below 
and in the next Section. We did not include a surface 
dependent magnetic depolarization contribution to 6iso, 
which depends on the magnetic susceptibility4^i^ 

Calculations were done using the DFT hybrid 
B3LYP)ii as well as the generalized gradient approx- 
imation (GGA), using the PW91 form3 The B3LYP 
calculations were found to yield generally better agree- 
ment with experiment fSec. lIIip than GGA. Site-centered 
GTO basis functions were associated with all the QM 
atoms. All-electron treatments were used for the O and 
Ti atoms, while the other QM atoms were represented 
using scalar-relativistic small core (RSC) pseudopoten- 
tials [also called effective core potentials (ECPs)]. All 
basis sets and ECPs listed below were taken from the 
EMSL website,^ O-ccntered (A4B20i5)i''" QM clusters 
were used in the calculations, suitably embedded with 
point charges and ePSPs, as discussed in Sec. Ill Bp . Test 
calculations for larger clusters were used to check size 
convergence, as discussed at the end of this section. 

The RSC ECPs for the QM atoms were Sr(28), Zr(28), 
Ba(46), Pb (60), where the number of core electrons is 
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shown in parenthesis. These pseudopotentials are gener- 
ally specified by the same label as their associated ba- 
sis sets listed below, except where otherwise indicated. 
For ePSPs fSec. Ill iTS)) . which have no associated GTOs, 
scalar-relativistic large core (RLC) ECPs were used as 
follows: CRENBS RLC for Ti(18) and Zr(36); Stuttgart 
RLC for Sr(36) and Ba(54). [The Ba ePSP was also 
used for the Pb ePSP atoms, because we could not 
find a large core Pb^^ quantum chemistry type pseu- 
dopotential in the standard databases.] The following 
GTO basis sets were found to give well-converged results: 
O(IGLO-III), Ti(cc-pwCVTZ-NR), Zr(cc-pwCVTZ-PP), 
and Pb(cc-pVTZ-PP). For the Sr and Ba atoms, the as- 
sociated Stuttgart RSC 1997 basis sets were used. The 
correlation consistent basis sets were employed in this 
study, because this basis function scries facilitates a sys- 
tematic study of chemical shift convergence with respect 
to basis set size. (We could not find correlation consistent 
basis sets for Sr and Ba.) The O(IGLO-III) basis set was 
used, because it was specifically designed for magnetic 
property calculations. In convergence tests, we found 
that it had quadruple-(^ (QZ) accuracy even though it 
has fewer basis functions. 

In SrTiOs, for example, B3LYP CSGT principle value 
components calculated with the aforementioned basis 
sets differed by no more than 10 ppm from results cal- 
culated with the O(cc-pwCVQZ), Ti (cc-pwCVQZ-NR), 
and Sr (def2-QZVP) basis sets, while Siso differed by only 
~ 8 ppm. In PbTiOs, changing the Pb basis set from DZ 
to TZ changed Siso by less than 10 ppm, while the differ- 
ence between inequivalent O sites changed by less than 
4 ppm. These results are consistent with calculations 
for SrTiOa, using two different gauge methods, shown in 
Table m where, in the infinite basis set limit, the CSGT 
and GIAO values should agree. The individual B3LYP 
components, ax = Uy and cr^, are seen to differ by 3 and 
26 ppm, respectively, while CTiso differs by only 6 ppm. 

Size convergence errors of about 30 ppm in the abso- 
lute value of the chemical shieldings were estimated from 
calculations on larger clusters, using the above basis sets. 
In SrTiOa, for example, increasing the 0-centercd clus- 
ter size from the Sr4Ti20i5 (21 atoms) to Sr4Tiio047 (61 
atoms) changed the RHF, B3LYP, and PW91 ctiso by -13, 
-27, and -25 ppm, respectively, while individual principal 
values changed by less than 16, 30, and 37 ppm. These 
test calculations used large core Ti pseudopotentials for 
the additional Ti atoms. 



III. RESULTS 

Calculated chemical shielding results for the prototyp- 
ical perovskites SrTiOs, BaTiOs, PbTiOs and PbZrOs 
are presented and compared to ^^O NMR single crys- 
tal and powder spectra chemical shift measurements in 
this section. Calculations were performed for embedded 
clusters corresponding to the following structural param- 
eters. For cubic SrTiOs, the lattice parameter of Ref. [49l 



was used. BaTiOa cubic and tetragonal P4mm struc- 
tures were taken from Ref. [50l : rhombohedral R3m from 
Ref. The PbTiOa tetragonal P4mm structure was 
taken from Ref. [5^. For PbZrOa, which has a complicated 
Pbam unit cell containing eight formula units, experi- 
mental lattice parameters from neutron scattering mea- 
surements were used together with internal coordinates 
determined from first principles calculations. — All cal- 
culations were carried out using the embedded cluster 
approach described in Sec. |TT1 using the hybrid B3LYP 
exchange-correlation functional. In some cases RHF and 
GGA/PW91 results are reported for comparison. 

Table U presents our calculated results for the principal 
values of the symmetrized chemical shielding tensor a. 
(The asymmetry in the chemical shielding tensors was 
less than 0.5 ppm; see Sec. Ill Al ) The principal axis frame 
of the target O atom is indicated by cc, y, z, where the z- 
axis is always identified with the local B-O-B direction. 
This is exact in the cubic and tetragonal crystals. For 
lower 0-site symmetry, the 2-axis well approximates the 
quasilinear B-O-B bond in all the structures considered 
here. 

A striking feature in Table|T]is the large ax — ay vs. a^ 
anisotropy of the principal values for each case, with two 
large negative (deshicldcd) principal values and one con- 
siderably smaller positive (shielded) cr^ principal value. 
The O atoms have their highest site-symmetry in the 
cubic crystals, resulting in an exact two-fold ax = ay 
degeneracy. For lower symmetry cases, the degeneracy 
is lifted, but the splitting remains small in most cases; 
the largest splitting is 65 ppm for the 01-4g site in 
PbZrOs. As shown by the GIAO calculations for SrTiOa, 
the anisotropy arises from the paramagnetic components. 
We note that the GIAO diamagnetic components in 
SrTiOa differ by only ~ 40 — 50 ppm from that of the the 
isolated, purely diamagnetic, 0^~ atom. 

In the fake SrTiOa calculation, where d-like {£ = 2) Ti- 
ccntered GTO basis functions were deliberately excluded 
from the calculation, the large x, y vs. z anisotropy is 
absent. The fake diamagnetic values are closer to the 
atomic 0^~ shielding, and the paramagnetic contribu- 
tions are close to isotropic. The x, y vs. z anisotropy of 
the calculated principal values and its relation to 0(2p) 
hybridization with the B-atom d-states is analyzed in 
Sec HVH 

Table [III compares the results for SrTiOa and BaTiOs 
in Table [T] with recent single crystal measurements^ 
and with earlier powder spectrum measurements of the 
isotropic shifts. Isotropic and uniaxial components of the 
chemical shift tensors are given, where (5ax = ((^z — <5iso)/2, 
using Eq. In SrTiOa, the calculated (5iso and Jax differ 
from experiment by ~ 20 ppm and ~ 9 ppm respectively. 
In BaTiOs, 5iso calculated and experiment values differ 
by at most 16 ppm. 

Table Iml compares the resuhs for PbTiOs and PbZrOs 
isotropic shifts ^iso (derived from Table |T| with well re- 
solved measured powder spectrai^ For both systems, the 
calculated results differ from experiment by no more than 
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21 ppm, with an average discrepancy of approximately 10 
ppm. Differences of relative shifts between chemical sites 
are smaller in most cases. For example, in tetragonal 
PbTiOs the calculated B3LYP splitting between the two 
incquivalcnt O sites is ~ 200 ppm, in good agreement 
with measured powder spectra.— 



IV. DISCUSSION 

A key feature in the electronic structure of transi- 
tion metal perovskites are the covalent interactions be- 
tween the 0(2p) and the (formally) unoccupied transition 
metal-d states. Indeed, the delicate balance between co- 
valent and electrostatic ionic interactions is responsible 
for the wide variety of interesting properties exhibited 
by these materials and related alloys. In this section, we 
analyze the calculated chemical shieldings in Sec. IIIII in 
relation to p-d hybridization. 

The measured NMR ^^O spectra show narrow well sep- 
arated peaks in these materials, indicating that second- 
order quadrupolar broadening and, thus, the electric field 
gradients (EFGs) are small^^ consistent with first princi- 
ples calculations of O EFGsj ^^i^^ For the systems studied 
here, all the O sites are clearly resolved in the measure- 
ments, and the chemical shielding tensor largely deter- 
mines the second order quadrupolar peak positionsi^ 



A. p-d hybridization and anisotropy of oxygen 
chemical shielding 

As mentioned, a striking feature in Table U is the 
large anisotropy of the chemical shielding principal val- 
ues. This is due to hybridization between the 0(2p) and 
virtual B-site d-states. The qualitative features can be 
understood from a simplified picture, which focuses on 
the B-O-B quasilinear structural unit (Fig. [T|). In a lin- 
ear molecule, the dependence of a on the direction of 
the applied field Bcxt is particularly simple, if we choose 
to locate the vector potential gauge origin at the target 
nuclcuSi ^^i^^ In this case, when Boxt is parallel to the 
molecular axis, only the diamagnetic (shielding) compo- 
nent contributes to ct; when Bext is perpendicular to the 
axis, both diamagnetic and paramagnetic (deshielding) 
components contribute. (Note that the calculated a are 
invariant under a change of the gauge origin. Conse- 
quently, we are free to interpret the u values as if they had 
been calculated with any particular choice of origin.) The 
diamagnetic component depends only on ground state 
wave functions, while the paramagnetic component de- 
pends on virtual transitions to unoccupied states, i.e., 
occupation of virtual states in the first order perturbed 
wave functions. In a DFT or HF calculation, this implies 
that large paramagnetic contributions could be expected, 
if there are low lying unoccupied one particle eigenstates 
that are strongly coupled to occupied states. This is the 
case in these materials, where there is strong coupling 



between 0(2p) and virtual B-site d-states. 

In cubic SrTiOa, for example, when Bext is applied 
along the Ti-0 bond direction (z-direction), the O nu- 
cleus is shielded by the applied field [i.e., the principal 
value (Tz = 90 is positive [Eq. ([T])]). By contrast, when 
Boxt is perpendicular to the Ti-0 bond, the O nucleus 
is strongly deshielded {a^ = cr^ = —343) as shown in 
Table HI According to this B-O-B picture, we would 
infer that only diamagnetic contributions should con- 
tribute to — 90; and the positive = 90 value is 
consistent with this. Similarly, we would attribute the 
deshielded = Oy = —343 value to paramagnetic con- 
tributions arising from the 0(2p)-Ti(3d) hybridization 
mechanism. This interpretation is supported by the fake 
SrTiOa calculations, where removal of Ti d basis func- 
tions quenches 0(2p)-Ti(3d) hybridization and results in 
nearly isotropic principal value components, <7x = ay ~ 

While qualitatively correct, the simple B-O-B model 
leaves out important crystalline effects. Thus the SrTiOs 
GIAO total shielding value Uz = 116 arises from cancel- 
lation of large diamagnetic a^^z — 371 and smaller para- 
magnetic ap^z = —255 components. Note that the fake 
ap^z = —262 value is similar to the normal crp,^ = —255 
value; in addition Cp^xy — <^-p,z in the fake calculation. 
We therefore attribute the presence of a non-zero ap^z to 
interactions with other atoms in the crystal, which are 
neglected in the simple B-O-B picture. 

Variations in B-O-B bond distances will affect 0(2p)- 
B(nd) hybridization, and this should be reflected by cor- 
responding changes in the calculated oxygen chemical 
shieldings. This is the case, as seen in Figure [21 which 
shows that the isotropic i5iso and uniaxial (Sax chemical 
shifts [Eq.([21)] exhibit a remarkably linear variation as 
a function of rg, the shortest B-0 bond length of the 
targeted O atomi^ Indeed, Siso changes by more than a 
factor of two over the plotted range. 

The linear behavior is largely due to the deshielding 
and Uy components in Table IH To frame the discussion 
in terms of the paramagnetic components of the shield- 
ing tensors, we first subtract the diamagnetic component 
of the isolated O^" atom, (t{0'^~) = 410 ppm, defining 
CTp,2 = - (t(0^~) and ap^xy = {(^x + o-y)/2 - cr(0^~). 
These are plotted in Fig. [31 An alternative choice would 
be to plot the GIAO paramagnetic component. Even had 
we done GIAO calculations for all the different 0-sites, 
this choice is not necessarily better, since both the dia- 
magnetic and paramagnetic GIAO components include 
contributions from neighboring atoms4^ Moreover, the 
separation into diamagnetic and paramagnetic compo- 
nents is, at best, only qualitatively useful, since the sep- 
arate components depend on the choice of gauge method. 
The present choice is physically motivated, using a well 
defined diamagnetic reference system. As indicated by 
the SrTiOs GIAO diamagnetic components, which differ 
by only c± 40 — 50 ppm from the O^- value (Tablc[I|, this 
subtraction largely removes the closed shell diamagnetic 
response. In any case, our definition of Up simply shifts 
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all (T principal values by a constant. As seen in Fig. [3l 
both ap^z and Cp^a^y vary linearly with rg. The average 
slope of (Jp^xy is ~ 5 times larger in magnitude than that 
of (Jp^z- This shows that the linear behavior in Fig. [2] is 
largely due to the variations of the deshielding ax and 
ay components, which arise from the p-d hybridization 
mechanism. 

The p-d hybridization mechanism also plays a key role 
in producing the anomalously large dynamical (Born) 
effective charge tensors Z*, which are universally seen 
in perovskite ferroelectrics for the O, transition metal 
B, and Pb atomsi^i^i^iSiH In cubic BaTiOg, SrTiOg, 
PbTiOa and PbZrOs, ^^^(O) (the component along the 
Ti-0 bond) takes on the values —5.59, —5.66, —5.83, 
and —4.81—, respectively, which is much larger than the 
nominal -2 value for O^". By contrast, the perpendic- 
ular components Z*^{0)=Z*y{0) are given by —2.11, 
—2.00; —2.56, —2.48, respectively, which arc much closer 
to the nominal value, since these stretch the B-0 bond 
only in second order. Similarly, all Z*(B'*+) ~ 6 — 7 
are anomalously large, since these involve the flow of dy- 
namic charge as the B-0 bond is stretched or compressed. 
[In PbTiOs and PbZrOa, Z*{Ph) ~ 4 is anomalously 
large compared to Z*(Ba) ~ Z*(Sr) ~ 2.6, due to strong 
Pb-0 covalcncy, and this is reflected in the somewhat 
more anomalous Z*^{0) and Z*y{0) in the Pb based 
crystals.] Artificially decreasing the p-d hybridization, 
as in Ref . [6^ resulted in nominal values of all the Z* . 

Similarly, the fake SrTiOa calculation eliminates B(d)- 
0(2p) paramagnetic contributions to a and largely re- 
moves the anisotropy between ax — ay and az com- 
ponents. B(d)-0(2p) covalency is essential for ferro- 
electricity in transition metal perovskites, mitigating 
the repulsion between otherwise rigid ion-cores, which 
tends to suppress ferroelectric distortion; under pres- 
sure, this mechanism eventually fails and ferroelectricity 
disappearsi^ Figure [3] shows that the the anisotropy be- 
tween ap^xy and ap^z increases linearly with decreasing r^, 
consistent with this picture, and indicating a strengthen- 
ing of the p-d hybridization mechanism as is reduced. 



B. Exchange and correlation effects 

Despite the fact that the B3LYP band gap is nearly 
a factor of two larger than that of GGA, the results in 
Tables U and [TTl are fairly similar, with the exception of 
PbTiOa, where the B3LYP results are in better agree- 
ment with experiment (Table IIIIl Sec. IIV C|) . Although 
larger band gaps result in larger energy denominators in 
the paramagnetic perturbative equations^^i^ there is lit- 
tle difference in the shieldings between the two methods. 
Similarly, a significant deshielded SrTiOs RHF paramag- 
netic component is evident in Tabled despite the much 
larger RHF HOMO-LUMO gap of 8.8 eV, compared to 
the B3LYP 3.6 eV gap. The apparent dependence of 
the chemical shielding on the band gap can be some- 
what misleading, however, since the transverse paramag- 



netic component of the first order current density could 
be made to vanish under suitable gauge transformations 
within the CSGT gauge method4^iS The entire shield- 
ing would then be given by the diamagnetic component, 
which depends only on the occupied single particle states. 
Similar observations were made regarding calculations of 
the the spontaneous polarization P and Z* in ferroelec- 
tric perovskites, where RHF results for these quantities 
were found to be in good numerical agreement with both 
experiment and with DFT calculations ^2*^ both meth- 
ods yielding anomalously large Z* . 



C. Comparison of oxygen chemical shifts with 
experiment 

The calculated results in Table |TT] for SrTiOs and 
BaTiOa are in very good agreement with single crys- 
tal measurements. Overall, both B3LYP and GGA yield 
similar agreement with experiment, while RHF is sig- 
nificantly worse. Comparisons for the cubic and ferro- 
electric P4mm tetragonal BaTiOa phases should keep 
in mind evidence for static and/or motional disorder 
of local structures suggested by experiment^i^SiZL and 
theory ; '^^i'^^ since the present calculations were performed 
for ordered crystalline structures. 

In PbTiOs, the Oax atom has two nn Ti atoms along 
the polar direction, with one short 1.77 A and one long 
2.39 A bond length, while the Ocq has two equal bonds 
of 1.98 A. The experimental assignment of the two mea- 
sured peaks is simplified by the fact that the relative 
integrated intensities of the two O spectral peaks corre- 
sponds to the 1:2 ratio of the Oeq to Oax atoms in the 
simple 5-atom primitive unit celli^ The large ~ 200 ppm 
experimentally observed^ splitting between the Oax and 
Ocq reflects their very different Ti-0 bond lengths. The 
B3LYP calculated chemical shifts in Table Hill reproduce 
the ~ 200 ppm splitting, while GGA underestimates it 
by 50 ppm. Basis set and cluster size convergence er- 
rors (Section III C|) can be expected to largely cancel in 
these calculated splittings, reducing the residual uncer- 
tainty to less than about 10 ppm. We are not aware 
of any measurements of the ^^O uniaxial asymmetry in 
PbTiOg. 

In PbZrOs the peak assignments are more difficult, 
because four of the five inequivalent oxygen sites have 
the same ratio of occurrence in the unit cell. In Ref. [5^ 
the peak assignment in the NMR spectrum was made 
based on the assumptions that (1) the proximity of the 
Zr cation plays the most important shielding role and (2) 
that the largest chemical shift corresponds to the largest 
bond length. The present calculations indicate that the 
first assumption is correct, but not the second. Instead, 
as seen in Fig. [21 the largest chemical shift corresponds 
to shortest bond distance between the targeted O atom 
and its nearest B atom. Indeed, the figure shows that 
this holds true, not only for PbZrOs, but across several 
compounds, over a wide range of chemical shifts. Our 
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calculations indicate that chemical shift site assignments 
for the 0(4) and 0(5) atoms in Ref. [H (their Table 1) 
should be reversed. Finally, we note that differences be- 
tween theory and experiment in the relative splittings are 
smaller than those of the absolute splittings, i.e., a rigid 
shift of the calculated B3LYP (GGA) values in Table Im] 
by ~ 10 ppm (~ —15 ppm) removes most of the the 
discrepancies. 

V. SUMMARY 

We have shown that first-principles embedded clus- 
ter calculations, using the DFT hybrid B3LYP method, 
can accurately calculate O chemical shifts for prototyp- 
ical perovskite structure transition metal oxides. Calcu- 
lated isotropic and uniaxial chemical shifts were found 
to be in good agreement with recent single crystal NMR 
^^O measurements for SrTiOa and BaTiOa. For ferro- 
electric P4mm BaTiOa and PbTiOa, the calculations 
accurately reproduced measured power spectra NMR 
isotropic chemical shifts Siso for the two inequivalent 
O sites. The large ~ 200 ppm experimental splitting 
in P4mm PbTiOs is well reproduced by B3LYP, but 
DFT/GGA underestimates ft by ~ 50 ppm. In PbZrOs, 
experimental and calculated (Sjso are in very good agree- 
ment, but experimental peak assignments are more diffi- 
cult, since four of the five inequivalent O sites appear in 
the same ratio. Our calculations, indicate a correction of 
the experimental assignments in Ref. IsBI . 

Our most notable findings are 1) a large anisotropy in 



the chemical shielding tensor, between deshielded ax — 
Uy and shielded Uz components, the latter principal axis 
being along the Ti-0 bond, and 2) a nearly linear vari- 
ation, across all the systems studied, of 5\so and (Sax as 
a function of B-O-B bond asymmetry. We have shown 
that the anisotropy and linear variation arise from large 
paramagnetic contributions to Ux and Uy due to virtual 
transitions between 0(2p) and unoccupied B(nd) states. 
A qualitative explanation of the anisotropy was given and 
then confirmed by calculations for a fake material with 
no d-states. 

We have shown that O NMR chemical shifts are a sen- 
sitive indicator of the local structure in perovskites with 
transition metal B-site atoms, due to covalent 0(2p)- 
B(nd) interactions. This indicates that ^^O NMR spec- 
troscopy, coupled with first principles calculations, can 
be an especially useful tool to study the local structure 
in complex perovskite alloys. 
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TABLE I: Calculated oxygen chemical shielding principal val- 
ues (in ppm) for SrTiOs, BaTiOa, PbTiOa and PbZrOg. In 
all cases, "z" denotes the principal axis that is most closely 
along the B-O bond direction. Unless otherwise indicated, 
the CSGT gauge method was used, with B3LYP for exchange 
and correlation. For some cases, GGA results are given in 
parenthesis. RHF and GIAO gauge method B3LYP results 
are given for SrTiOs. GIAO results for the isolated diamag- 
netic 0^~ atom are also shown. See the text for a discussion 
of the fake SrTiOs calculation. For PbZrOs, we adopt the 
site labeling convention of Ref. [s^ the 0(i) site notation of 
Ref. [5^ is also given in brackets. 









ay 


O'z 






SrTiOa (cubic) 








-343 (-353) 


-343 (-353) 


on /' 

yu ( 'ioj 






347 


347 


o t ± 


GIAO paramag. 




-694 


-694 


-255 


GIAU total 




-347 


-347 


116 


RHF 




-134 


-134 


201 






410 


410 


410 


fake ca-lculation 














397 


397 


448 


v^i/^vy paiaiiiag. 




-228 


-228 


-262 


LriAU total 




169 


169 


186 






BaTiOs 






cubic 




-413 (-414) 


-413 (-414) 


97 (49) 


P4mm 





ax -505 (-480) 


-505 (-480) 


116 (70) 




0, 


,q -401 (-407) 


-353 (-355) 


87 (40) 


R3m 




-423 (-416) 


-412 (-406) 


97 (51) 






PbTi03(P4mm) 









ax -608 (-562) 


-608 (-562) 163 (123) 




0, 


-263 (-286) 


-219 (-228) 


13 (-32) 






PbZr03(Pbam) 




01-4g [0(1)] 




-172 (-197) 


-107 (-130) 


85 ( 48) 


01'-4g [0(2)] 




-132 (-158) 


-104 (-126) 


77 ( 40) 


02-81 [0(3)] 




-147 (-171) 


-139 (-162) 


100 ( 64) 


03-4f [0(4)]'' 




-93 (-118) 


-64 ( -90) 


62 ( 26) 


04-4e [0(5)] 




-251 (-269) 


-223 (-239) 175 (136) 



"Note: our calculations suggest that the assignment of the 03- 
4f[0(4)] and 04-4e[0(5)] peaks should bo reversed in Ref. [ssl. as 
done here. 
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TABLE II: Comparison between calculated and measured single crystal and powder spectra ^^O chemical shifts (in ppm) for 
BaTiOs and SrTiOa. Calculated isotropic shifts (5iso and uniaxial components 5ax are derived from Table HI referenced to liquid 
water, o-i^o*""^ = 287.5 ppm (Ref. [t^. Calculated values are from B3LYP (GGA values are in parenthesis), and RHF results are 
also shown for SrTiOs. 







Theory 


Expt. 


<5ax 

Theory 


Expt. 


SrTiOa 












cubic 




486 (507) 


467 ± 5^^, 465^ 


-144 (-132) 


-135.3 ± 5^* 


cubic RHF 




310 




-112 




BaTiOs 

cubic 
P4mm 


Oax 
Oeq 


530 (547) 
585 (584) 
510 (528) 


546 ± 5° 
570 ± 5", 553^ 564'' 
520 ± 5°, 530^ 523'' 


-170 (-154) 
-207 (-183) 
-155 (-140) 


-150 ± 1" 
-171 ± 1" 
-142 ± 1" 



"Single crystal experimental values are from Ref. 1 54 
''Powder results Ref. 
'^Powder results Ref. 
■^Powder results Ref. TTI 



TABLE III: Comparison between calculated and experimen- 
tal isotropic chemical shifts Siso (in ppm) for PbTiOa and 
PbZrOa. Calculated isotropic shifts are derived from Table U 
referenced to liquid water, aZT'' = 287.5 ppm (Refill). Cal- 
culated values are from B3LYP (GGA values are in parenthe- 
sis). For PbZrOa, we adopt the site labeling convention of 
Ref. [5^: the 0(i) site notation of Ref. Issl is also given. 



Theory 



Expt.'' 



PbTiOa 

Oax 
Oca 



639 (621) 
444 (469) 



647 ± 2 
443 ± 2 



PbZrOa 

01- 4g [0(1)] 
01'-4g [0(2)] 

02- 8i [0(3)] 

03- 4f [0(4)] 

04- 4e [0(5)] 



352 (380) 
340 (369) 
349 (377) 
319 (348) 
387 (412) 



365 ± 2 
351 ± 2 
356 ± 2 
329 ± 2 
408 ± 2 



"Experimental values are from Ref. [55| 

''Note: our calculations suggest that the assignment of the 03- 
4f[0(4)] and 04-4e[0(5)] peaks should be reversed in Ref. Issl. as 
done here. 
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FIG. 1: Illustration of cluster embedding for the O-centered A4B2O15 quantum atom (QM) cluster (Sec. IIIB Ifl . shown for a 
[110] plane of the ideal perovskite structure. The QM A, B, and O atoms are depicted by filled black, blue and red circles, 
respectively. In the left panel, the QM atoms are embedded in point charges, "+" or "— " signs colored coded according to the 
QM atom they replace in the crystalline lattice. In the right panel, the boundary O QM atoms have had their nearest and 
next-nearest neighbor cation point charges replaced by "empty" pseudopotentials (Sec. Ill iTS)) . indicated by thick large circles, 
with corresponding color coding. 
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FIG. 2: Calculated oxygen isotropic Siso and uniaxial Jax chemical shifts in SrTiOa (ST), BaTiOa (BT), PbTiOa (PT), PbZrOa 
(PZ), as a function of Vs, the shortest B-0 bond length of the targeted O atom. The notation B--0--B indicates O atoms 
with two equidistant nn B atoms, and B-0--B indicates an O atom with one short and one long nn B bond {e.g., along the 
polar c-axis in P4mm PbTiOa). The straight lines are linear fits. 
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FIG. 3; Calculated paramagnetic (Tp,z and (J^,xy components (see text), as a function of r^, the shortest B-0 bond length of 
the targeted O atom. The straight lines are linear regressions of the points. Symbols as in Fig. [2] 



